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The Eddy Dissipation Rate is an important metric for representing the intensity of 
atmospheric turbulence and is used as an input parameter for predicting the decay of 
aircraft wake vortices. In this study, the forecasts of eddy dissipation rates obtained from 
the current state-of-the-art mesoscale model are evaluated for terminal area applications. 
The Weather Research and Forecast mesoscale model is used to simulate the planetary 
boundary layer at high horizontal and vertical mesh resolutions. The Bougeault-Lacarrere 
and the Mellor-Yamada-Janjic schemes implemented in the Weather Research and Forecast 
model are evaluated against data collected during the National Aeronautics and Space 
Administration’s Memphis Wake Vortex Field Experiment. Comparisons with other 
observations are included as well. 


I. Introduction 

D ECAY and transport of aircraft wake vortices depend on the strength of the ambient atmospheric turbulence, 
among other factors. Current fast-time wake transport and decay prediction models (Sarpkaya et al. 2001; 
Holzapfel and Robins 2004; Proctor et al. 2006; Proctor 2009) use eddy dissipation rate (EDR) as the parameter that 
controls the decay of wake vortices due to atmospheric turbulence. The EDR can be obtained from wind data 
measured at high frequency with various sensors such as sonic anemometers. However, the cost of deploying a 
comprehensive array of sensors that can adequately cover the entire spatial extant of an airport terminal area is 
prohibitive and therefore measurements can be taken at the most at few points in space. Another disadvantage is the 
problem of providing forecasts of EDR for several hours following the measurement time. Weather forecasts for 
supporting the terminal area operations can be obtained with the help of numerical models and substantial progress 
has been made in mesoscale modeling over the last three decades. The operational models include the National 
Center for Atmospheric Research’s (NCAR) Weather Research and Forecasting (WRF) model (Klemp et al. 2007), 
Science Applications International Corporation’s Operational Multiscale Environment model with Grid Adaptivity 
(Bacon et al. 2000), the United States Navy’s Coupled Ocean/ Atmosphere Mesoscale Prediction System (Hodur 
1997), the Mesoscale Atmospheric Simulation System (MASS) model (Meso, Inc. 1995), the Mesoscale Model 5 
(MM5) developed jointly by the Pennsylvania State University and the National Center for Atmospheric Research 
(Dudhia 1993) and the National Centers for Environmental Prediction’s (NCEP) Rapid Update Cycle (RUC) 
developed by Benjamin et al. (2004). 

There are many challenges in using mesoscale models for wake vortex applications. First, the computational 
cost of running the models at an appropriate spatial and temporal resolution for wake applications can be 
prohibitive. A fine mesh-resolution can be achieved in the vertical, but a practical resolution in the horizontal is 
limited to between 4 km and 10 km. The restriction on horizontal mesh resolution is due to computational constraints 
as well as the validity of various physics parameterizations used in the model. This horizontal resolution can, at 
best, resolve the whole terminal area by just a few grid points (Figure 1). The use of coarse grid resolution for wake 
applications can introduce uncertainties in the forecast of winds and turbulence quantities. Secondly, the physics 
parameterizations in mesoscale models are scale specific which makes it difficult to simulate complex atmospheric 
phenomenon that are inherently multiscale. The grid resolution related uncertainties are in addition to the 
uncertainties due to the nonlinearity of the governing equation set solved by the models. Another uncertainty in 
model forecasts is due to inadequate initial and boundary conditions. 
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Some recent applications of mesoscale modeling to support terminal area operations in general and wake vortex 
research in particular include: Kaplan et al. (1999), Chamey et al. (2000), Kaplan et al. (2006), Freeh et al. (2007), 
and Ringley et al. (2007). These studies have had mixed results in terms of success in predicting winds, temperature 
and turbulence quantities. Temperature predictions usually had the least amount of error, whereas, the EDR 
obtained from the turbulence kinetic energy (TKE) predictions were either significantly overestimated or 
underestimated. 

In this study, the National Center for Atmospheric Research (NCAR) state-of-the-art mesoscale model is 
evaluated as a tool for providing forecasts of EDR. The model forecasts are compared against data observed at 
Memphis International Airport during the month of August in 1995. This data was collected during the Memphis 
Wake Field Experiment which was sponsored by the National Aeronautics and Space Administration (NASA) 
during the month of August in 1995. In the following sections the Memphis Field Experiment is described briefly 
and the results of the WRF simulations are discussed in detail. 



Figure 1. Memphis International Airport is shown in the left panel and its approximate length and width is 
marked on the figure. The horizontal mesh resolution of mesoscale models is usually on the order of a few 
kilometers and at best can represent the entire airport by a few points. In the vertical, however a high- 
resolution can be prescribed. The grid shown in the right panel has a resolution of 40 m in the boundary 
layer. Above the boundary layer the grid is stretched all the way to the top of the computational domain. 

II. Eddy Dissipation Rate 

According to Kolomogrov’s hypothesis, the large eddies in the atmosphere are inherently unstable. These large 
eddies break up and transfer their energy to smaller-sized eddies. This transfer or “cascade” of energy continues 
from larger to smaller eddies, until viscous dissipation becomes dominant and the TKE is converted into heat. The 
rate at which the TKE is dissipated is known as the EDR. In other words, EDR is the rate at which the TKE is 
absorbed by breaking down eddies into smaller and smaller eddies until it is converted into heat by viscous forces 
(Figure 2). The length scale of the large energy containing eddies, / 0> also known as the integral scale of the 
turbulence, can be derived as follows: 


3/2 


l 0 CC- 


where, s is the eddy dissipation rate in nr/s 3 . The turbulence kinetic energy, e in nr/s 1 is given by: 


2/„2 . 




e = — \u' 2 +v' 2 +w' 2 


( 1 ) 


( 2 ) 
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In Eq. (2), u\ v\ w' are the turbulent fluctuations of u , v, and w-velocity components from some time-averaged 
mean value and the bar symbol denotes time averages. 

Re-arranging Eq. (1) gives: 


s = 


Ce 3/2 


(3) 


where, C is a proportionality constant of order one. The eddy dissipation rates can be estimated directly from high 
frequency wind measurements (e.g., Campbell et al. 1997) and from similarity theory (Han et al. 2000). 

In mesoscale models, EDR can be diagnosed from the model’s TKE closure schemes that are used for computing 
subgrid scale processes. The accuracy of the TKE predictions and EDR calculations depend on the underlying 
assumptions and approximations in the particular TKE closure scheme used by the model. Since the TKE is a 
standard output of most mesoscale models, the unknowns left in Eq. (3) are the integral length scale, / 0 , and the 
proportionality constant, C. Eq. (3) is sometimes further simplified as follows: 

3/2 

^ = — (4) 


where the proportionality constant, C, and the integral length scale, / 0 , have been combined into a length scale, Z e , 
that needs to be specified. 

Gerz et al. (2005) have developed the Nowcasting Wake Vortex Impact Variables (NOWVIV) model to forecast 
the cross winds and turbulence quantities for supporting wake vortex studies. The NOWVIV tool is built upon the 
MM5 model with a level-3 turbulence closure scheme (Burk and Thompson 1989). In one of their studies, Freeh et 
al. (2007) have evaluated the EDR predictions obtained from the NOWVIV model. The inner-most nest in their 
study had a horizontal grid resolution of 2.1 km and a stretched vertical grid that placed 26 of the 60 vertical levels 
below 1100m. Freeh et al. (2007) defined L e equal to a constant 311m and found that resulting EDR values were 
underestimated above an altitude of 100m. They suggested using direct model outputs of EDR for heights above 
100m and TKE-derived estimates below 100m. Ringley et al. (2007) assumed L e = 336m (instead of 31 lm), which 
was based on their analysis of 40m wind observations. In their study, the EDR was overestimated by one or two 
orders of magnitude near the surface, whereas the predictions at an altitude of 40m compared well with observations. 
Kaplan et al. (2006) used the MASS model to run simulations on a fine-scale mesh resolution that assumed a 
horizontal grid spacing of 167m and 60 vertical levels below 2 km. In Kaplan et al. (2006), the length scale, T e , was 
calculated as a function of Richardson number (Therry and Lacarrere 1983; and Bougeault and Lacarrere 1989) and 
the EDR estimates compared well with the observations. Charney et al. (2000) used similarity theory (Han et al. 
2000) to estimate EDR and TKE. This methodology, however, significantly overestimated the EDR. 


Transfer of energy to smaller scales 



Figure 2. The Kolomogrov Energy Spectrum. Energy range for eddies of different sizes and their associated 
length scales are shown in the figure as a function of wavelength. 
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III. The Memphis Wake Vortex Field Experiment 

A comprehensive field experiment to measure wake vortices and the ambient meteorological conditions was 
conducted at the Memphis International Airport in Memphis, Tennessee from August 6 through August 29, 1995 
(Campbell, et al. 1997). The experiment was sponsored under NASA Langley Research Center’s Aircraft Vortex 
Spacing System (AVOSS) project (Hinton 1995; Perry et al. 1997). Several sensor systems were deployed to collect 
detailed meteorological data to study the effect of atmospheric conditions on the behavior of wake vortices. The 
sensors included radiosondes, sodars, a wind profiler, one 150 ft high meteorological tower, a radio acoustic 
sounding system (RASS) and NASA Langley’s OV-10 research aircraft. The radiosondes were used to measure 
winds and temperature measurements (10s averages) at 50m vertical resolution. The OV-10 aircraft was flown at 
selected times and took measurements of temperature and winds at a sample rate of 10 Hz. Temperature (5 min 
averages) was measured using RASS every 30 min at 97 vertical levels from 145m to 1492m. The 150 ft (45.7m) 
meteorological tower (located at latitude 35.029167°N and longitude 89.981 111°W and an elevation of 97.37m 
above sea level) was equipped with a large array of sensor systems (Table 1). Winds, temperature and moisture 
were measured from the tower at 5m, 10m, 20m, 30 m, and 42m heights. Turbulence quantities (TKE and EDR) 
were estimated from wind measurements at 5m and 40m heights. Rain rate, soil temperature, soil moisture, 
barometric pressure and incoming and outgoing solar radiation also were measured by the sensors deployed on the 
meteorological tower. In addition to the tower data, standard meteorological data such as atmospheric pressure, 
temperature, moisture, cloud cover, and visibility, etc. was obtained from the National Weather Service’s Surface 
Aerodrome Observations (SAO) and the Automated Surface Observations System (ASOS). Figure 3 shows some of 
the sensors deployed at the Memphis International Airport during the field experiment. 



Figure 3. Memphis Wake Vortex Field Experiment. A B757 landing in the background of the NASA lidar 
van (top left); profiler (top right); Michael Mathews (left) and Dr. Fred Proctor (right) standing next to the 
radiometer mounted at the base of the meteorological tower (bottom left); and Dr. Tim Dasey in the lidar van 
(bottom right). 
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Table 1: Data Collection by the 150 ft Meteorological Tower 


Measured Quantity 

Height (m) 

/(Hz) 

Averaging Time 

Sensor System 

Accuracy 

temperature (°C) 

5, 10, 20,30, 42 

1 

1 min 

SAVPAK 

0.3°C 

relative humidity (%) 

5, 10, 20,30, 42 

1 

1 min 

SAVPAK 

3% 

wind speed (m/s) 

5, 10, 20,30, 42 

1 

1 min 

SAVPAK 

2% 

wind direction (degrees) 

5, 10, 20,30, 42 

1 

1 min 

SAVPAK 

1 degree 

virtual temperature (K) 

5,40 

10 

1 min 

FLUXPAK 

0.05°C 

w-velocity (m/s) 

5,40 

10 

1 min 

FLUXPAK 

0.05m/s 

v-velocity (m/s) 

5,40 

10 

1 min 

FLUXPAK 

0.05m/s 

w-velocity (m/s) 

5,40 

10 

1 min 

FLUXPAK 

0.05m/s 

pressure (mb) 

2 

1 

1 min 

BAROMETER 

0.1 5mb 

net radiation (W/m 2 ) 

2 

1 

1 min 

RADIOMETER 

n/a 

incoming radiation (W/m 2 ) 

2 

1 

1 min 

RADIOMETER 

n/a 

outgoing radiation (W/m 2 ) 

2 

1 

1 min 

RADIOMETER 

n/a 

soil temperature (°C) 

1 0cm below ground 

1 

1 min 

SOILPAK 

0.05°C 

soil moisture (%) 

10cm below ground 

1 

1 min 

SOILPAK 

n/a 

accumulated rain (inches) 

2 

n/a 

1 min 

SOILPAK 

O.Olin 


In order to compare the model predictions with the observed data, a period was selected for the numerical 
simulations from August 15 through August 20, 1995, which coincided with the times that observed data was 
collected during AVOSS. During the August 1995 field experiment, the weather conditions followed the typical 
summertime pattern for the Memphis area. The days were hazy, hot, and humid with generally light winds 
(Campbell et al. 1997; Zak 1995). The atmosphere in the nighttime was hazy, warm and very humid with light 
winds (Campbell et al. 1997; Zak 1995). The observed daily maximum and minimum temperatures, dew point 
temperatures, and maximum wind speed at 5 m were obtained from SAVPAK and are listed in Table 2. The sunrise 
and sunset times for the period of interest are also listed in Table 2. The values in Table 2 reflect this pattern with 
the exception of August 19. The large wind speed recorded on August 19 is due to a gust front that moved through 
the airport. Precipitation was recorded on August 17 and August 19, from thunderstorms that developed nearby 
during the late afternoon and the evening hours. 


Table 2: Some Observations at Memphis, Tennessee 


date 

sunrise 

sunset 

min T(°F) 

max T(°F) 

min T d (°F) 

max T d (°F) 

max wspd (m/s) 

19950815 

0619CDT 

1949CDT 

76.19 

92.12 

69.78 

76.78 

4.75 

19950816 

0620CDT 

1948CDT 

76.46 

94.37 

71.71 

77.12 

4.67 

19950817 

0621 CDT 

1947CDT 

77.99 

96.06 

72.95 

80.15 

6.89 

19950818 

0621 CDT 

1945CDT 

77.63 

97.14 

71.04 

77.59 

6.33 

19950819 

0622CDT 

1944CDT 

73.36 

95.79 

66.68 

80.45 

18.38 

19950820 

0623 CDT 

1943 CDT 

SAVPAK data was not available for this date 


Central Daylight Time (CDT) = Coordinated Universal Time (UTC) - 5 hr 
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IV. Mesoscale Modeling 

The WRF model (Klemp et al. 2007) is a mesoscale weather prediction model that has been developed by NCAR 
with help from various partners in academia and federal research laboratories. It has governing equations for the 
time-dependent, three-dimensional, nonhydrostatic, fully compressible Navier-Stokes equations. The model grid 
structure is based on Arakawa C-grid staggering (Arakawa and Lamb 1977). In the vertical, WRF uses a terrain- 
following pressure formulation, where the top of the domain is a constant pressure surface. The mesh can be 
stretched in the vertical to provide higher resolution in the boundary layer. One-way and two-way nesting options 
are available to provide high horizontal grid resolution in regions of interest. The horizontal and vertical advection 
terms can be discretized anywhere from second to sixth-order spatial accuracy. Fifth-order upwind differencing in 
the horizontal direction and third-order upwind differencing in the vertical direction are recommended by the 
developers of the WRF model. The time integration is done using an explicit third-order Runge-Kutta time 
marching scheme, with smaller time steps for acoustic and gravity wave modes. Monotonic as well as positive- 
definite schemes can be used for the advection of microphysical scalars, chemical scalars and TKE. High resolution 
terrain and other land use datasets are part of the standard WRF distribution. 

Extensive physics packages for modeling the surface layer interactions, boundary layer turbulence processes and 
cumulus parameterizations have been implemented in the model. Twelve different cloud microphysics schemes of 
varying complexity are available, as well as parameterizations for radiation physics (four different schemes for 
longwave radiation and seven schemes for shortwave radiation). The user has the choice of six schemes for surface 
layer interactions. Nine different schemes have been implemented in the WRF model for the parameterization of 
sub-grid scale diffusion processes in the planetary boundary layer (PBL). The schemes vary in complexity and 
range from simple non-local-K schemes (e.g., Troen and Mahrt 1986) to large eddy simulation (LES) for very high 
resolution grids. The LES option in WRF has not been extensively tested or used. There are a number of TKE- 
based closures implemented in the WRF model (Janjic 1994; Bougeault and Lacarrere 1989; Nakanishi and Niino 
2004). The TKE-closure schemes implemented in many mesoscale models are usually some variant of the Mellor- 
Yamada scheme (Mellor and Yamada 1982), in which the prognostic equation for the TKE is computed. 

In this study, the Mellor- Yamada- Janjic scheme (Janjic 1994) was used, which is a variant of the Level 2.5 
Mellor-Yamada turbulence scheme. The Mellor-Yamada scheme in its original form is singular in the limit of 
vanishing shear. The Mellor- Yamada- Janjic (MYJ) scheme is an improved version of the Mellor-Yamada scheme 
that remedies the singularity problem. In addition to the MYJ PBL scheme, the Bougeault and Lacarrere (BL) 
scheme was also evaluated. The BL scheme is 1.5 order and solves the prognostic equation for the TKE. The 
specification of the length scale and other constants in the BL PBL scheme is discussed in Bougeault and Lacarrere 
(1989). The MYJ PBL scheme is described in detail in Janjic (2001). When these two PBL schemes are used, the 
WRF model output consists of TKE and the master length scale which are then used to compute the EDR in the 
post-processing step. 


V. Simulations Results 

In this section the results of WRF simulations are discussed in detail. Six days of simulations were performed 
from August 15 to August 20, 1995 using the two different PBL schemes. Prior to these simulations, a detailed 
sensitivity study was conducted to evaluate the effect of various WRF input parameters. The results were compared 
with available surface and upper air observations. 

A. Simulations Setup and Summary 

Three levels of grid nesting were utilized in the simulations in order to achieve acceptable grid resolution while 
maintaining computational efficiency. The WRF simulation domain consisted of an outer most domain bounded 
between 124.06°W and 55.93°W in the longitude and 18.08°N and 49.14°N in the latitude with a horizontal mesh 
resolution of 3 6km. Two higher resolution nests were defined within the outermost domain. Domain 2 was 
bounded between 111.54°W and 68.45°W in the longitude and 24.84°N and 43.86°N in the latitude with a 
horizontal mesh resolution of 12 km. The innermost domain was bounded between 99.74°W and 79.97°W in the 
longitude and 30.66°N and 39.29°N in the latitude with a horizontal mesh resolution of 4 km. High-resolution (1 km) 
terrain and land use datasets were used for the innermost domain while lower resolution terrain was used in the 
coarser nests. The WRF computational domain and the terrain for the innermost domain are shown in Figure 4. 

The simulations were initialized using North American Regional Reanalysis (Mesinger et al. 2006) data from the 
National Climate Data Center (NCDC). The fifth-order upwind-biased scheme was used in the horizontal and the 


6 

American Institute of Aeronautics and Astronautics 



third-order upwind-biased scheme was used in the vertical within a three-stage Runge-Kutta explicit time-marching 
scheme. The monotonic scheme option was used for the transport of both the microphysical scalars and the TKE. 

Several sensitivity runs were made for August 16, 1995. These runs were conducted to evaluate the effects of 
spin-up time, time step, time step for radiation calls, vertical mesh resolution, PBL scheme, and the microphysics 
scheme. Two different spin-up times were tested in the simulations: 1200UTC and 1800UTC on August 15, 1995. 
The simulations were run for 36 and 30 hours depending on initialization time for a 24 hr forecast. The two different 
spin-up times of 6hr and 12 hr were found to make little difference and therefore, a spin-up time of 6hr was used. 

The time step based on the CFL criteria can be set according to: 


At 

ai max ^ [— 

V 3 U max 


( 5 ) 


where, C cf i is the CFL criteria. If the fifth-order upwind scheme is used within a three-stage explicit Runge-Kutta 
time marching scheme then this constant is set to 1.42. The wind speed magnitudes in the jet stream can be 
substantial and therefore, u max is usually set to 100mA to insure stability. The developers of WRF also suggest using 
the time step according to the following relation: 


At ma x= c i Ax ( 6 ) 

where, c x = 6s/km is a constant, Ax is the grid resolution in kilometers and At is in seconds (e.g., a Ax of 18 km would 
suggest a At = 108s). In practice, however, these formulations may not work every time and the time step may need 
to be further reduced depending on the resolution of the vertical grid, underlying terrain, etc. In the current study it 
was found that the time step chosen by using both Eq. (5) and Eq. (6) was not sufficiently small enough. In some 
cases the simulation became unstable and stopped without giving any errors. The optimal time step required to 
maintain stability had to be determined by trial and error. In the simulations described in this section a time-step of 
15s was used. 

The Rapid Radiative Transfer Model (RRTM) longwave radiation scheme (Mlawer et al. 1997) and the Goddard 
shortwave scheme (Chou and Suarez 1999) were used to parameterize the effects of both longwave and shortwave 
radiative transfer in the atmosphere. The radiation time step was set to 4 min (the suggested value for the radiation 
time step is 1 min per 1 km of grid size, e.g., it should be 10mm if Ax = 10 km). Different radiation time step values 
were tested and it was found that a higher frequency call to both longwave and shortwave radiation schemes 
produced a better simulation of the diurnal cycle. 

The schemes tested for microphysical processes were the WRF Single-Moment 5 -class scheme (details are given 
in the WRF User’s Guide), Lin et al. (1983) scheme, and the Thompson et al. (2004) scheme. The Thompson 
microphysics scheme performed better compared to the Lin microphysics scheme. There were instances when the 
Lin microphysics scheme generated precipitation over Memphis when no precipitation was observed on ground. 
The use of WRF Single-Moment 5-class scheme also degraded the accuracy of the forecasts. The cumulus 
parameterization schemes were switched off for the simulations. 

The BL and MYJ PBL schemes were tested using uniform vertical mesh resolutions of 40m, 20m, and 10m 
below 800m. The vertical mesh was stretched above 800m to coarser resolutions. The vertical grid levels that 
corresponded to these mesh resolutions were 60, 90, and 125, respectively. Increasing the vertical mesh resolution 
did not necessarily result in improved results. In some cases an increase in the vertical resolution degraded the 
quality of overall predictions. 

The surface time histories of temperature and dew point temperature were noisy when the MYJ PBL scheme was 
used. Increasing the value of the 6 th order filter coefficient suppressed the noise to a certain degree, but not 
completely. It should be noted that the higher-order schemes (both centered and upwind-biased) implemented in the 
ARW WRF model are neither TVD nor positive-definite. The transport of microphysical scalars and TKE however 
can be either positive-definite or monotonic. 

The forecasts obtained from the sensitivity studies were compared with observations and their accuracy was 
quantified in terms of root mean square error, mean absolute error, and bias. The input configuration that yielded a 
sufficiently accurate forecast for Memphis for August 16, 1995 was used for all other days (Appendix A). Two sets 
of simulations were conducted (from August 15 to August 20) with all input parameters held constant except the 
PBL scheme. The results of these simulations are discussed in the next two sections. 


7 

American Institute of Aeronautics and Astronautics 




110°W 100°W 90°W 80°W 70°W 



Topography height (meters MSL) 


[ in 1 i i J i 


10 160 310 460 610 760 910 1060 1210 1360 


Figure 4. The top panel shows the WRF computational domain with higher resolution nests and the bottom 
panel shows the terrain contours for the innermost 4 km resolution nest. 
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B. Comparison of WRF Temperature and Winds with Observations 

The simulation results were compared with surface observations and upper air soundings at three locations: 
Memphis, Tennessee, Nashville, Tennessee and Little Rock, Arkansas. The station information is given in Table 3 
and Table 4. The accuracy of simulations was quantified in terms of root mean square error (Error n ns ), mean 
absolute error (Error m ae ) and Bias: 


Em™ 



(7) 


Errors 


n 



obs 


( 8 ) 


Bias = -Y j (xf -x° bs ) 


(9) 


The sensitivity runs were conducted for August 16, 1995, since mesoscale models are known to have difficulty 
in accurately simulating strongly stable nocturnal boundary layers with weak synoptic forcing (Hanna et al. 2010), 
as in this case. 

The model preformed reasonably well and was able to simulate the main features of the diurnal cycle for August 
16, 1995. The transition of the nighttime stable boundary layer to the daytime convective boundary layer was 
effectively simulated as well as the formation of the lower level nocturnal jet (Figure 5). With some exceptions, the 
surface temperature and wind predictions from the WRF model were in good agreement with observations at 
Memphis for the different dates. 



time (utc) time (utc) 


Figure 5. WRF time history at Memphis for August 16’ 1995. The left panel shows the time history of 
temperature (°F) and the right panel shows the wind speed (m/s). MYJ PBL Scheme was used. 


The WRF simulation results using the MYJ PBL scheme are compared with the surface observations from the 
Memphis airport and from the meteorological tower (SAVPAK) in Figures 6-7. The errors in temperature and wind 
velocities computed for the two sets of simulations at three different locations are tabulated in Tables 5-10. The 
errors at Memphis for wind velocities were fairly low, regardless of the choice for PBL schemes. The errors in 
temperature and dew point temperature were generally larger at Nashville and Little Rock for both PBL schemes 
compared to errors at Memphis. The errors in wind velocities at Nashville and Little Rock were generally low and 
similar to the errors at Memphis. 
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A closer look at the simulations suggests that the larger forecast errors present in some of the simulations (e.g., 
August 18) can be attributed to inadequate initialization data. These errors can be reduced by using better data 
assimilation and initialization techniques (e.g., observation nudging) and or a different value of spin-up time. In the 
case of August 18 simulation, the initialization was not improved by increasing the spin-up time to 12 hr, as seen by 
the dashed line in Figure 6 representing 12 hr spin-up time. 

In some cases, the model failed to capture local subgrid scale features such as convective activity that resulted in 
large forecast errors. On August 17, for example, the simulated temperature and dew point temperature compared 
well initially with the observations, but failed to capture the sharp temperature drop at approximately 2300UTC 
(Figure 6). The 2m airport sensor and the 5m SAVPAK sensor show a drop of more than 10 degrees, which is not 
forecasted by the model. 

This discrepancy between model forecast and observations was further investigated. The surface observations at 
Memphis recorded a thunderstorm between 2215UTC and 2323UTC with rainfall starting at 2215UTC and ending 
at 2310UTC. Maximum wind gusts of 9.2 m/s (18te) were recorded at the surface. The NEXRAD data over the 
Memphis airport at 2136UTC showed radar reflectivity (RRF) values as high as 60 dBZ. Figure 8 shows the 
composite of lower level RRF values at different times over the Memphis airport. The echo tops were as high as 
13.7 km (45 kft) in the vicinity of Memphis during this period. Similarly, on August 19, a thunderstorm developed 
between 2031UTC and 2141UTC with rainfall starting at 2107UTC and ending at 2132UTC. Maximum wind gusts 
of 22.1mA (43 kts) were recorded at the surface. The model was again not able to simulate this convective event. 
The sharp temperature drop and the increase in the observed surface winds due to the moving gust front were not 
captured in the model simulation (Figures 6-7). Further studies can be done to evaluate other microphysics schemes 
implemented in the WRF model. An increase in horizontal mesh resolution may help in resolving the small scale 
convective features. 


Table 3: Surface Observation Stations 


Location 

Station ID 

Latitude 

Longitude 

Elevation 

Memphis, TN 

MEM 

35.042°N 

89.976°W 

103.9m 

Nashville, TN 

BNA 

36.124°N 

86.678°W 

182.6m 

Little Rock, AR 

LIT 

34.729°N 

92.224°W 

79.9m 


Table 4: Upper Air Soundings 


Location 

Station ID 

Latitude 

Longitude 

Elevation 

Memphis, TN 

MEM 

35.029°N 

89.981°W 

97.37m 

Nashville, TN 

BNA 

36.25°N 

86.55°W 

210m 

Little Rock, AR 

LZK 

34.72°N 

92.23°W 

78m 


The comparisons of upper air soundings from WRF and Memphis on August 16, 1995 are shown in Figures 9- 
12. Also shown are the comparison of the predicted WRF wind speed and direction with the profile generated by 
Lincoln Labs. The Lincoln Labs profile, in which data measured from several different types of instruments was 
fused together to generate a single profile, agrees well with the model simulations (Figures 9-10). The low-ascent 
field experiment sounding data which was at a much higher resolution is plotted for Memphis in Figure 1 1 . 
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Figure 6. Comparison of WRF simulations with observations. Observed temperatures at 2m for Memphis 
are shown as squares and the SAVPAK 5m data is given by small diamonds in the figure. The SAVPAK data 
frequency is every one minute. WRF simulation results are denoted by solid lines. The SAVPAK data was 
not available for August 20, 1995. 
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Figure 7. Comparison of WRF simulations with observations. Observed winds at 10 m for Memphis are 
depicted by squares and the 5m SAVPAK data is given by small diamonds in the figure. The SAVPAK data 
frequency is every one minute. WRF simulation results are denoted by solid lines. The SAVPAK data was 
not available for August 20, 1995. 
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Table 5: Memphis, Tennessee (Bougeault-Lacarrere Scheme) 


date 

T(°F) 

T„(°F) 

u(m/s) 

v(m/s) 

rmse 

mae 

bias 

rmse 

mae 

bias 

rmse 

mae 

bias 

rmse 

mae 

bias 

950815 

1.17 

0.96 

- 0.27 

2.08 

1.69 

0.77 

1.25 

1.04 

- 0.72 

0.54 

0.45 

- 0.20 

950816 

1.53 

1.25 

0.96 

1.28 

0.99 

-0.58 

1.52 

1.19 

0.52 

0.92 

0.71 

-0.31 

950817 

3.05 

1.74 

1.62 

1.60 

1.25 

- 0.99 

2.03 

1.70 

1.22 

1.64 

1.31 

- 1.26 

950818 

3.80 

3.07 

3.07 

1.15 

0.96 

- 0.14 

1.36 

1.12 

1.07 

1.61 

1.24 

- 1.00 

950819 

4.03 

2.19 

2.17 

2.20 

1.87 

- 1.14 

1.60 

1.21 

- 0.47 

1.65 

1.31 

0.07 

950820 

4.76 

4.38 

3.54 

3.06 

2.71 

1.69 

1.97 

1.63 

- 0.04 

2.78 

2.35 

0.55 


August 16, 1995 (950816) was used for the different sensitivity runs. 


Table 6: Nashville, Tennessee (Bougeault-Lacarrere Scheme) 


date 

T(°F) 

T d (°F) 

u(m/s) 

v(m/s) 

rmse 

mae 

bias 

rmse 

mae 

bias 

rmse 

mae 

bias 

rmse 

mae 

bias 

950815 

2.52 

2.24 

1.86 

2.86 

2.53 

- 2.31 

1.75 

1.56 

1.11 

2.05 

1.81 

0.63 

950816 

2.43 

1.96 

1.88 

2.28 

1.82 

1.35 

2.03 

1.58 

- 0.31 

1.86 

1.55 

0.57 

950817 

1.84 

1.55 

1.49 

1.07 

0.83 

0.59 

1.94 

1.51 

- 0.80 

1.41 

1.15 

0.79 

950818 

5.95 

3.58 

3.26 

2.49 

2.22 

2.17 

2.05 

1.65 

- 1.42 

2.52 

1.66 

0.30 

950819 

6.34 

5.80 

4.59 

3.68 

3.29 

3.29 

2.34 

1.90 

- 1.65 

2.30 

1.95 

- 0.37 

950820 

2.81 

2.47 

2.47 

4.84 

4.12 

3.39 

2.59 

1.87 

- 0.79 

2.55 

2.06 

0.05 


Table 7: Little Rock, Arkansas (Bougeault-Lacarrere Scheme) 


date 

T(°F) 

T„(°F) 

u(m/s) 

v(m/s) 

rmse 

mae 

bias 

rmse 

mae 

bias 

rmse 

mae 

bias 

rmse 

mae 

bias 

950815 

2.39 

2.07 

- 0.85 

4.10 

3.08 

2.36 

1.90 

1.61 

- 0.62 

0.79 

0.60 

- 0.01 

950816 

2.03 

1.61 

- 0.91 

2.76 

2.44 

0.73 

2.04 

1.82 

0.54 

1.29 

1.05 

- 0.08 

950817 

2.89 

2.25 

- 1.43 

3.62 

2.78 

1.48 

2.75 

2.31 

2.03 

2.11 

1.81 

- 0.03 

950818 

2.57 

1.61 

- 1.35 

1.92 

1.69 

- 0.68 

2.84 

2.70 

2.11 

1.97 

1.57 

- 0.49 

950819 

4.50 

3.10 

- 3.07 

2.52 

1.84 

0.56 

3.11 

2.60 

1.33 

2.07 

1.69 

- 0.06 

950820 

3.81 

3.38 

2.14 

4.27 

3.90 

3.83 

1.44 

1.21 

0.11 

2.29 

1.72 

0.91 
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Table 8: Memphis, Tennessee (Mellor-Yamada-Janjic Scheme) 


date 

T(°F) 

T„(°F) 

u(m/s) 

v(m/s) 

rmse 

mae 

bias 

rmse 

mae 

bias 

rmse 

mae 

bias 

rmse 

mae 

bias 

950815 

1.43 

1.23 

- 0.76 

1.42 

1.02 

0.11 

1.41 

1.21 

- 0.76 

0.58 

0.48 

0.21 

950816 

1.40 

1.15 

0.81 

1.37 

1.02 

-0.66 

1.68 

1.32 

0.74 

0.88 

0.74 

-0.43 

950817 

2.56 

1.43 

1.11 

1.98 

1.63 

- 1.46 

2.15 

1.84 

1.43 

1.68 

1.36 

- 1.33 

950818 

3.87 

3.07 

3.07 

1.38 

1.04 

- 0.84 

1.25 

0.98 

0.86 

1.78 

1.40 

- 1.12 

950819 

4.92 

2.80 

2.79 

2.76 

2.34 

- 2.30 

1.73 

1.12 

- 0.48 

1.60 

1.21 

- 0.04 

950820 

4.39 

3.53 

3.46 

3.43 

3.05 

1.41 

2.48 

2.03 

- 0.16 

4.26 

3.09 

1.43 


August 16, 1995 (950816) was used for the different sensitivity runs. 


Table 9: Nashville, Tennessee (Mellor-Yamada-Janjic Scheme) 


date 

T(°F) 

T„(°F) 

u(m/s) 

v(m/s) 

rmse 

mae 

bias 

rmse 

mae 

bias 

rmse 

mae 

bias 

rmse 

mae 

bias 

950815 

2.24 

1.94 

1.58 

3.24 

2.87 

- 2.85 

1.82 

1.66 

1.43 

2.01 

1.75 

0.63 

950816 

2.22 

1.83 

1.55 

2.26 

1.91 

0.59 

1.95 

1.44 

0.01 

2.17 

1.83 

1.16 

950817 

1.89 

1.60 

1.54 

0.84 

0.69 

- 0.28 

1.88 

1.47 

- 0.91 

1.65 

1.35 

0.79 

950818 

5.78 

3.43 

3.08 

1.50 

1.29 

0.93 

2.03 

1.64 

- 1.38 

2.61 

1.90 

0.21 

950819 

6.85 

6.31 

3.22 

3.45 

3.25 

3.25 

3.48 

2.47 

- 2.34 

2.42 

2.20 

- 0.38 

950820 

3.13 

2.64 

2.02 

5.42 

4.57 

4.32 

2.25 

1.78 

- 0.98 

2.05 

1.73 

0.04 


Table 10: Little Rock, Arkansas (Mellor-Yamada-Janjic Scheme) 


date 

T(°F) 

T d (°F) 

u(m/s) 

v(m/s) 

rmse 

mae 

bias 

rmse 

mae 

bias 

rmse 

mae 

bias 

rmse 

mae 

bias 

950815 

2.43 

2.00 

- 1.04 

3.70 

2.92 

2.04 

1.70 

1.42 

- 0.36 

0.92 

0.70 

0.002 

950816 

2.25 

1.75 

- 1.13 

2.21 

1.93 

0.55 

2.10 

1.78 

0.84 

1.23 

0.99 

- 0.04 

950817 

3.04 

2.35 

- 1.93 

3.24 

2.42 

1.21 

2.91 

2.49 

1.94 

2.18 

1.84 

- 0.01 

950818 

2.61 

1.96 

- 0.80 

1.47 

1.29 

- 0.05 

2.90 

2.63 

2.37 

2.03 

1.65 

- 0.69 

950819 

5.18 

3.60 

- 3.46 

1.91 

1.42 

1.27 

2.98 

2.45 

2.01 

1.88 

1.38 

0.007 

950820 

3.17 

2.73 

0.63 

4.03 

3.66 

3.54 

1.86 

1.51 

- 0.09 

3.00 

2.35 

1.05 


14 

American Institute of Aeronautics and Astronautics 





Figure 8. Memphis NEXRAD data for August 17, 1995. NEXRAD Site ID: KNQA. The composite of lower 
level radar reflectivity at different times is shown in the figure. The start of convective activity over the 
Memphis International Airport is after 2119UTC. Areas of high radar reflectivity (> 60 dBZ) can be seen in 
the vicinity of the airport. 
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Figure 9. Comparison of WRF simulation data with observations (August 16, 1995). Wind speed profile for 
Memphis at time = 6 hr into the simulation is shown in the left panel and at time = 9 hr into the simulation is 
shown in the right panel. Points are observations and solid lines are WRF computations. Red points in the 
right panel correspond to the profile generated by Lincoln Labs which fused the data measured from several 
different types of instruments. 
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Figure 10. Comparison of WRF simulation data with observations (August 16, 1995). Wind direction 
profile for Memphis at time = 6hr into the simulation is shown in the left panel and at time = 9 hr into the 
simulation is shown in the right panel. Points are observations and solid lines are WRF computations. Red 
points in the right panel correspond to the profile generated by Lincoln Labs which fused the data measured 
from several different types of instruments. 
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temperature (F) temperature (F) 

Figure 11. Comparison of WRF simulation data with observations (August 16, 1995). Vertical 
temperature profile for Memphis at time = 6hr into the simulation is shown in the left panel and at time = 9hr 
into the simulation is shown in the right panel. Points are observations and solid lines are WRF 
computations. 




temperature (F) temperature (F) 

Figure 12. Comparison of WRF simulation data with observations (August 16, 1995). Vertical 
temperature profile for Nashville at time = 6 hr into the simulation is shown in the left panel and for Little 
Rock at time = 6 hr into the simulation is shown in the right panel. Points are observations and solid lines are 
WRF computations. 
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C. Comparison of Model derived EDR with measured EDR 

The EDR is diagnosed from WRF model output of TKE and the master length scale. Previous studies have 
shown that the use of model predicted master length scale in EDR calculation results in large errors. Some 
researchers have suggested using a constant value of length scale for EDR calculation, therefore, in addition to the 
EDR calculated using the model output of the master length scale, three different values of a constant length scale 
were also tested. The constant values used were 336m as suggested by Ringley et al. (2007); 311m suggested by 
Freeh et al. (2007) and 300m. The use of different values of a constant length scale did not change the EDR values 
significantly. The errors in the TKE and EDR computed for the two sets of WRF simulations (using the BL and the 
MYJ PBL schemes) are given in Tables 11-12. The EDR in the tables was calculated using both the constant length 
scale (300m) as well as the model output of length scale. In general the errors were smaller when the MYJ PBL 
scheme was used and errors were smaller when a constant value of the length scale was used instead of the model 
output of the length scale. The errors listed in Tables 11-12 were computed using only the data points between 
1200UTC and 2400UTC which roughly corresponds to the daytime well-mixed convective boundary layer. The 
FLUXPAK data was not available for the daytime hours of August 15, 1995. On August 17 and August 19, 1995, 
the spike in both TKE and the EDR occurs around the time of convective activity which was not captured by the 
model. This resulted in large errors in TKE and EDR predictions on these two dates for both PBL schemes. 

Compared to the RMS errors in EDR predictions reported by Freeh et al. (2007) which were on the order of 
0.004m 2 A 3 , the RMS errors in EDR using MYJ scheme for the current study are an order of magnitude lower when 
using a constant length scale. Ringley at al. (2007) also conducted an extensive evaluation of the MASS model 
predictions of TKE and EDR against the 1997 AVOSS data collected at the Dallas-Fort Worth Airport. They did 
not report quantitative measures such as RMSE but instead provided plots of model predicted TKE and EDR 
compared with measurements. The plots given in Ringley et al. (2007) indicate a model underestimation of TKE by 
as much as 75% in the nighttime boundary layer. There are also large errors in EDR predictions during both the 
daytime and the nighttime boundary layer. In their evaluation of mesoscale models for dispersion applications, 
Hanna et al. (2010) noted that the TKE and other boundary layer parameters were “least well simulated on very 
stable nights” when the models default to the threshold values of TKE. The Nonhydrostatic Mesoscale Model 
(WRF-NMM) and the MM5 model were used in their study. For the daytime boundary layer, the TKE predictions 
were usually within a factor of two of the observations with degradation in nighttime predictions. On average they 
found that the WRF-NMM and MM5 predictions underestimated the TKE by 20-50%. 

The WRF simulation results of TKE and EDR using the MYJ PBL scheme are compared with the met tower 
observations (FLUXPAK) in Figures 13-14. The WRF prediction for the daytime convective boundary layer 
matches the observed data fairly well. However, the results for the nighttime stable boundary layer were 
overestimated (solid black lines in Figures 13-14). An increase in the vertical mesh resolution resulted in an 
improvement for the nighttime prediction of EDR by the BL PBL scheme. In the case of MYJ PBL scheme, 
however, an increase in the vertical mesh resolution slightly improved the results for the nighttime predictions but 
degraded the quality of the overall forecast. In general, grid convergence could not be achieved for either scheme. 

The simulation results showed that during the nighttime stable boundary layer the model often defaulted to a 
threshold value for the TKE which was too large. The threshold values are typically used in the absence of 
sufficient instability and wind shear to ensure non-zero TKE. This resulted in large values of EDR during the 
nighttime stable boundary layer (solid black lines in Figure 14). For wake applications this would imply a 
prediction of rapid vortex decay. Janjic (personal communications) suggested using a lower value for the TKE 
threshold. The floor value of TKE was decreased in the WRF code from 10" m Is to 10 ' m Is and the simulations 
were performed again. The results for the modified code are depicted by the dashed orange lines in Figures 13-14. 
The results for August 16, 1995 showed substantial improvement in both the TKE and EDR predictions and for 
other days gave a conservative prediction for EDR. In general, the decrease in the TKE threshold resulted in either 
an approximately average value or an underestimation of EDR for the nighttime stable boundary layer. The 
conservative prediction of EDR in this case was based on a realistic TKE threshold value and is also desirable for 
operational safety concerns. A lower value of EDR implies larger wake decay times and adds toward the factor of 
safety in the fast-time wake model predictions. 

Finally it should be noted that although this study was conducted to evaluate the two PBL schemes in the WRF 
model, the strong coupling of various physics and dynamics modules with each other makes it difficult to isolate the 
performance of one specific module or scheme. The current study reflects the performance of the two PBL schemes 
as implemented in the ARW version of the WRF model. Further studies can be conducted to evaluate other 
mesoscale models such as the WRF-NMM from NCEP and the Operational Multiscale Environment model with 
Grid Adaptivity (OMEGA). 
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Table 11: Memphis, Tennessee (Bougeault-Lacarrere Scheme) 


date 

TKE(m 2 /s 2 ) 

EDR Lw (m 2 /s 3 ) 

EDR Lc (m 2 /s 3 ) 

rmse 

mae 

bias 

rmse 

mae 

bias 

rmse 

mae 

bias 

950816 

0.69 

0.56 

0.51 

0.005 

0.004 

0.004 

0.004 

0.003 

0.003 

950817 

1.16 

0.74 

0.22 

0.006 

0.004 

0.003 

0.005 

0.003 

0.002 

950818 

0.53 

0.40 

0.29 

0.005 

0.003 

0.003 

0.004 

0.003 

0.002 

950819 

1.71 

0.88 

-0.02 

0.02 

0.007 

-0.0005 

0.02 

0.006 

-0.001 

950820 

0.72 

0.56 

0.25 

0.007 

0.005 

0.005 

0.005 

0.004 

0.004 


EDR Lw = EDR calculated using the length scale from model output 
EDR Lc = EDR calculated using a constant length scale of 300 m 


Table 12: Memphis, Tennessee (Mellor-Yamada-Janjic Scheme) 


date 

TKE(m 2 /s 2 ) 

EDR Lw (m 2 /s 3 ) 

EDR Lc (m 2 /s 3 ) 

rmse 

mae 

bias 

rmse 

mae 

bias 

rmse 

mae 

bias 

950816 

0.35 

0.19 

-0.07 

0.0009 

0.0007 

-0.0007 

0.0006 

0.0004 

0.0003 

950817 

1.18 

0.52 

-0.50 

0.002 

0.001 

-0.001 

0.002 

0.001 

-0.0006 

950818 

0.36 

0.27 

-0.25 

0.001 

0.001 

0.0007 

0.0005 

0.0003 

-0.00003 

950819 

1.58 

0.54 

-0.42 

0.02 

0.005 

-0.002 

0.02 

0.004 

-0.003 

950820 

0.58 

0.40 

-0.32 

0.001 

0.001 

0.001 

0.001 

0.0008 

0.0006 
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Figure 13. Comparison of simulated TKE values with observed data. Simulation results using the Mellor- 
Yamada-Janjic (MYJ) PBL scheme are shown for different dates. Points represent observations taken within 
three different time windows (5 min, 15min and 30 min). WRF computations are given by lines. Az=40/w. 
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Figure 14. Comparison of simulated EDR values with observed data. Simulation results using the Mellor- 
Yamada-Janjic (MYJ) PBL scheme are shown for different dates. Points represent observations taken within 
three different time windows (5 min, 15min and 30 min). WRF computations are given by lines. Az=40/w. 
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VI. Summary 

Several simulations were conducted using different physics options and model run configurations. The WRF 
model was able to capture the main features of the complete diurnal cycle such as the development of a low level jet 
during the nighttime stable boundary layer and the transition of the nighttime stable boundary layer into the daytime 
convective boundary layer. The forecasts of surface temperature and winds compared favorably both qualitatively 
and quantitatively with observations. The cases with large errors in the model forecasts were either due to poor 
initialization or due to model’s inability to capture local subgrid-scale convective activity. 

The TKE and the EDR predictions compared well with observations for the daytime convective boundary layer, 
but large errors were observed for the nighttime stable boundary layer. The BL PBL scheme underestimated the 
EDR during the nighttime and the MYJ PBL scheme overestimated EDR during the nighttime. An increase in the 
vertical mesh resolution did not give consistent improvement in predictions and in some cases degraded the 
accuracy of the overall forecast. In general, grid convergence could not be achieved for either scheme. Decreasing 
the value of TKE threshold for the MYJ PBL scheme resulted in more realistic EDR predictions for the nighttime 
boundary layer. Overall, the errors were smaller when the MYJ PBL scheme was used compared to the BL PBL 
scheme and errors were smaller when a constant value of the length scale was used instead of the model predicted 
length scale. 

The mesoscale models have proved invaluable for meteorological forecast and analysis and for detailed forensic 
reconstruction of events of interest. The errors in temperature and wind forecasts have been reduced significantly 
over the years but in their current state they are not recommended for real-time predictions of turbulence parameters 
such as EDR to support terminal area operations. Further research leading to improvements in PBL closure schemes 
is needed, especially for accurate simulations of the nighttime stable boundary layer with weak synoptic forcing. 

Appendix A 

The main input parameters used in the simulations are listed in this appendix. These parameters were held 
constant for the two sets of simulations using the Mellor-Yamada-Janjic (MYJ) and the Bougeault and Lacarrere 
(BL) schemes. The values used in some of the sensitivity studies for August 16, 1995 are shown in bold font. 

• Time-marching: 3 rd order Runge-Kutta 

• Time step: 15s 

o 20s, 30s, 60s 

• Momentum advection: 5 th order upwind in the horizontal and 3 rd order upwind in the vertical 

• Ax = Ay = 4 km 

• Number of vertical levels: 60 

o 90, 125 

• Vertical mesh resolution (A z) near surface: 40m 

o 20/w, 10 m 

• Longwave Radiation Scheme: RRTM 

• Shortwave Radiation Scheme: Goddard 

• Radiation call: 4 min 

o 2/wi7f, 10 min 

• Microphysics Scheme: Thompson 

o Lin, WRF Single-Moment 5-class 

• Advection of TKE and microphysical scalars: Monotonic 

o Positive-definite 

• Model initialization using the NARR data from NCDC 

• Spin-up Time: 6hr 

o 12 hr 
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